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Nonlinear response of quantum cascade structures 
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The gain spectrum of a terahertz quantum cascade laser is analyzed by a non equilibrium Green's functions 
approach. Higher harmonics of the response function were retrievable, providing a way to approach nonlinear 
phenomena in quantum cascade lasers theoretically. Gain is simulated under operation conditions and results 
are presented both for linear response and strong laser fields. An iterative way of reconstructing the field 
strength inside the laser cavity at lasing conditions is described using a measured value of the level of the losses 
of the studied system. Comparison with recent experimental data from time-domain-spectroscopy indicates 
that the experimental situation is beyond linear response. 



Possible coherent radiation in the terahertz range 
has been a very strong motivation for research in the 
field of Terahertz Quantum Cascade Lasers 1,2 (THz- 
QCLs), which would enable a wide range of applica- 
tions such as imaging 3 and spectroscopy 4 . However 
compact devices operating over cryogenic temperatures 
are a practical requirement for applications and cur- 
rently the most promising designs are based on resonant 
phonon extraction 5 , achieving operating temperatures up 
to ~ 200 K. 6 

The key physical quantity in any QCL is the gain which 
describes the amplification of the optical field in the het- 
crostructurc material. In recent years this quantity has 
been measured in detail in Time-Domain-Spectroscopy 
(TDS) experiments 7,8 where THz-QCLs are probed by 
ultra short pulses providing information on both phase 
and amplitude of the transmitted pulse, whereafter the 
gain spectrum is reconstructed by a Fourier transform. 
The pulse is made as strong as possible in order to get 
a good signal to noise ratio but it is not known how the 
system dynamics are affected by such a measurement. 
The simulation of THz QCLs relies on a consistent treat- 
ment of tunneling and scattering, either by hybrid den- 
sity matrix/rate equation schemes 9-13 or more evolved 
Non Equilibrium Green's Function (NEGF) theory 14-17 . 
Here we present an extension of our NEGF scheme 18 to- 
wards the treatment of high intensities inside the QCL, 
going beyond linear response to an external electromag- 
netic field. 

In this article we consider a time-dependent electric 
field F(t) = Fdc + F ac cos(ftt) in the cavity, which re- 
flects both the applied bias (fdc) and the electric compo- 
nent (F ac ) of a monochromatic field of angular frequency 

in the cavity. Going beyond linear response, this re- 
quires the solution of the time-dependent Kadanoff-Baym 
equation for the lesser and retarded Green's functions, 
G< j3 (k,t 1 ,t 2 ) and G I ^(k,t 1 ,t 2 ), respectively 19 . Here 
a, p denote the states in growth direction and k the in- 
plane momentum. The periodicity in time allows for a 
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Fourier decomposition of the Green's functions 
G{k;h,h) = y [dEj2 e- iE ^- t ^ h G h (k,E)e- i 

^ J h 

(1) 

and similarly for the self-energies. This provides a set of 
equations for the Green's functions for given self-energies 
£ Q ^.h(k, E), which arc defined analogously. This proce- 
dure follows essentially the concepts outlined in Ref. 20 
and details will be given elsewhere. Here the terms with 
h = correspond to the stationary transport considered 
before , while the higher order terms take into account 
the ac field. For the fields considered in this manuscript 
we used h = —2,— 1,0, 1,2, while checking that increas- 
ing \h\ did not change the results (generally higher val- 
ues of \h\ are required with increasing ratios eF^d/Ml, 
where d is the period of the QCL structure). Relations 
to observables are made through the h = and h = 1 
components although the higher orders effect the lower 
ones implicitly. The Green's functions allow for a deter- 
mination of the current, where G^a h (k,E) provides the 
dc current for h = and the ac current with frequency 
51 for h = 1. Dividing the latter by F ac provides the 
conductivity, directly related to the gain coefficient. 

Here we consider the sample studied in Ref. 8 as shown 
in Fig. 1(a). In Fig. 1(b) we show the current-bias rela- 
tion, calculated by our model, together with experimen- 
tal data. The original experimental data refer to the 
bias U along the entire heterostructure, containing 175 
periods of length d as well as contact regions, where ad- 
ditional bias drops. In a similar experiment 21 a voltage 
drop of around 3.0 V as well as a 517 series resistance 
were assumed for the data analysis. Here, we use a volt- 
age drop of 3.8 V in the contacts for converting the ex- 
perimental bias to eFa c d as displayed in Fig. 1(b) and 
Fig. 3(b). 

Fig. 1(b) displays simulations with different interface 
roughness scattering parameters in order to calibrate one 
of the parameters used, namely the average (RMS) of the 
roughness height 77. This enters in the matrix element for 
the interface roughness scattering self energy 22 together 
with the typical size A = 10 nm of the roughness layers. 
According to Fig. 1(b) it is clear that a lower value of 
77 suppresses the current flow, and regarded as a fit pa- 
rameter 77 = 0.20 nm is a better value. 23 All subsequent 
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FIG. 1. (Color online) (a) QCL structure considered with 
the main states contributing to its operation and (b) calcu- 
lated current-bias relation for r\ = 0.2 nm (solid) and rj = 0.3 
nm (dotted-dashed) together with experimental data of Ref. 8 
(dashed). The shift in current at simulated operating condi- 
tions (dotted) is discussed later. The marked points at differ- 
ent current densities are analyzed in Fig. 2. 



simulations where carried out with rj = 0.20 nm and a 
lattice temperature of 77 K entering the occupation of the 
phonon modes. The experimental heat-sink temperature 
was ~ 33 K 8 , but the lattice temperature for resonant 
phonon extraction THz-QCLs is typically higher 24 due 
to heating effects. 

We note that for both simulations the main peak as 
well as the low-bias behavior arc in good agreement with 
the experimental data. In contrast, for bias drops per 
period of about 40-45 mV, we observe a spurious extra 
peak due to tunneling between state 1 and state 3 of the 
next neighboring period. Such extra peaks for long-range 
tunneling have been observed in our model for other THz- 
structures as well, and we currently attribute them to the 
lack of electron-electron scattering processes 25 . 

Fig. 2 shows the calculated gain spectra (lines) for 
three different dc biases, as depicted in Fig. 1(b), both 
for low (a) and high (b) ac field strength. For compar- 
ison, wc display the experimental data from Ref. 8 in 
both parts. The low current point (□/--) is considered 
to be at the bias where state 1' and 4 become aligned 
and electrons start to tunnel through the system as those 
levels are at resonance. Here we are far below threshold 
and mainly absorption is seen at the laser frequency of 
2.2 THz, as state 4 - the lower laser state - is populated, 
while 5 is still mostly empty. This is changed as current 



FIG. 2. (Color online) Gain spectrum at different operation 
conditions for small eF ac d = 1 meV (a) and larger eF ac d = 
6 meV (b) ac field strength. Open symbols are experimental 
data, and lines are simulated data. Green (dotted-dashed) 
corresponds to simulations at 34 mV/period, red (dashed) 
to 50 mV/period and blue (solid) to 54 mV/period. The 
experimental currents were 63, 319 and 403 A/cm 2 for the 
low, medium and high current respectively. 

is increased and we approach the medium current point 
(o/--). This is taken where the system is almost at, but 
still below the threshold current. Here the states 1' and 
5 start to align, creating population inversion at the laser 
frequency. At operating conditions, where the third and 
last point (A/-) at high current is taken, gain is above 
the level of the losses 8 of 18 cm" 1 and can now sustain 
lasing as the population inversion is at its maximum. 

The simulations at low and high ac field strength 
shown in Fig. 2 differ drastically, but the general picture 
is that around the laser frequency, large losses at low 
current develop into gain as bias is raised, as observed 
by Refs. 8 and 26. At a more detailed level, the low ac 
field strength simulations exhibit stronger features which 
are not reflected in the experimental data. At higher 
ac field strength however, the overall agreement becomes 
better, mainly due to a redistribution of carriers (bleach- 
ing). However, the strong absorption feature around 1.2 
THz for medium and high current density in the experi- 
mental data does not appear in the simulations, for which 
we do not have any explanation. The better agreement 
of experimental data with high ac field, indicates that 
the experimental conditions are beyond linear response. 

Here it is important to address the fundamental dif- 
ferences between experiment and simulation. In the sim- 
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ulations, a monochromatic ac field is applied and the 
gain spectrum is constructed frequency by frequency. In 
the experimental case, the situation is quite different. A 
pulse, containing all frequencies within the bandwidth 
(typically 3 THz 7 ) is sent into the sample, and the way 
this pulse has changed by passing through the structure 
determines the gain spectrum. Compared to the simu- 
lations, where we only measure at the frequency where 
we excite, this is the opposite, as all frequencies are sub- 
ject to excitations and all frequencies are also measured. 
Therefore, our modeling can only be seen as approximate. 
In addition, the experimental data is the difference of 
measurements between the unbiased QCL and the QCL 
at the chosen measuring bias. This way the background 
is effectively subtracted. If the structure exhibits less 
losses at some point than it does at zero bias, this is 
measured as gain. In the simulations however, we only 
extract the gain from the conductivity extracted from the 
h = I component of the Green's functions, as we do not 
have to take losses into account and thus only look at 
the intrinsic gain spectra. Simulations at very low bias, 
i.e. the off-state, show absorption peaks at 0.9 THz and 
2.7 THz, which could explain corresponding features in 
the experimental gain spectra. 

We have demonstrated, that the response varies signif- 
icantly with the ac field strength. Thus it is important 
to question whether the ac field strengths used in the 
simulations are comparable to their experimental coun- 
terpart. In order for the effects of high ac field strengths 
shown in Fig. 2(b) to be of any relevance, the power 
coupled into the QCL structure during the experimental 
measurements must be sufficient. Addressing this ques- 
tion, consider the experimental situation governing the 
in-coupling of light 8 : A pumping femtosecond pulse of 
125 mW hits the emitter section of the same QCL as 
described in Fig. 1. The pulse generates a photocurrent 
giving an electric field transient that is coupled across 
an air distance of 4 /urn into the QCL section of interest. 
Our value eF ac d = 6 meV corresponds to a power of 40 
mW in the cavity, which requires an extremely efficient 
conversion in the emitter and good coupling between the 
structures. It is far from clear, whether the probing field 
can reach these intensities. Strong ac fields at lasing con- 
ditions would be capable of generating these effects, but 
this would then only contribute above threshold current. 

In a working laser, the gain will clamp at the level 
of the losses, as the population inversion will stabilize 
around the configuration where the inversion lost to the 
optical field will be balanced by the injector efficiency. 
This happens at different intensities for different bias 
points, and using this fact the intensity at gain clamping 
and thus the power in the QCL can be reconstructed the- 
oretically. One can use the level of the losses measured 
in Ref. 8 in order to iteratively calculate the intensity at 
operating conditions. Fig. 3(a) shows gain spectra at a 
bias of 54 mV/period that have been simulated at various 
ac-field strengths. It is clear from Fig. 3 that higher ac 
field strengths effectively lower both gain and absorption, 
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FIG. 3. (Color online) (a) Gain spectra at a bias of 54 
mV/period at different ac field strengths, (b) Calculated in- 
tensities in the waveguide (dashed) and corresponding exper- 
imental data (full line). 



and give rise to gain bleaching. For all bias points above 
threshold simulations were carried out at different ac field 
strengths in order to see which ones gave gain matching 
the level of the losses at 18 cm" 1 . The corresponding in- 
tensities were then extracted and are shown in Fig. 3(b). 
For the maximum intensity of 0.08 mW//im 2 the wave 
guide area of 800 /im 2 provides a power of 64 mW inside 
the QCL. This seems reasonable taking into account that 
only a part is coupled out through the mirror. 

To show the importance of including higher orders of 
the Fourier decomposed Green's function in Eq. (1), sim- 
ulations with \h\ < 1 only and also \h\ < 3 is shown in 
Fig. 3(a) for eF ac d = 8 meV which is the highest ac- 
field used. \h\ < 1 (dashed) shows a clear deviation from 
the \h\ < 2 case (full line) while the simulations with 
\h\ < 3 confirm the quality of our \h\ < 2 calculations for 
eF ac d < /Ml. 

The -F ac -assisted current at the intensities shown in 
Fig. 3(b) is displayed as a dotted line in Fig. 1(b). It 
increases proportional to the intensity compared to the 
non-lasing current. Thus when lasing sets in, we see a 
kink in the current, just as in the experimental data of 
Fig. 1(b) at 53 mV per period. As the calculated kink ap- 
pears somewhat stronger, our calculated lasing intensities 
could be a little bit too high. This may be related to the 
fact, that the experimental lasing frequency of 2.2 THz 
is not precisely at the peak of the gain spectrum. 

In conclusion, we have simulated gain under opera- 
tion by including higher orders of the Fourier decomposed 
Green's function in order to include nonlinear effects. We 
have found a way to calculate the power of the laser us- 
ing an experimental value of the level of the losses and 
by iteratively matching the gain to that level and then 
extract the intensity of such a configuration. It has also 
been shown that gain bleaches under high intensity con- 
ditions and that this might be a non negligible effect in 
THz-TDS measurements. 
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